function V = V(globvar,par,na1,rb,Vini,fmat,fprob)

s=par{1};beta=par{2};theta=par{3};ra=par{4};fb=par{5};
mar=par{6};infi=par{7};diffini=par{8};difftol=par{9};

amat  = globvar{1};
zmat  = globvar{2}; 
ztran = globvar{3};
pmat  = globvar{4};   
ptran = globvar{5};
bmat  = (1-theta)*exp(pmat);

% vectorize state space [b a z p f] (nh*nb*na*nz*nf)
nb = length(bmat);
na = length(amat);
nz = length(zmat);
np = length(pmat);
nf = length(fmat);
ns = nb*na*nz*np*nf; 

% statespace = [bvec avec zvec pvec];
bvec    = repmat(bmat,na*nz*np*nf,1);
avec    = repmat(repelem(amat,nb),nz*np*nf,1);
zvec    = repmat(repelem(zmat,nb*na),np*nf,1);
pvec    = repmat(repelem(pmat,nb*na*nz),nf,1); 
fvec    = repelem(fmat,nb*na*nz*np,1); 

% choice set of a1
a1choice  = repelem(linspace(0,max(amat),na1),ns,1);  % ns*na1, same for every state

% points to be interpolated
a1interp  = repmat(vec(a1choice'),1,nz*np*nf);

b1ninterp = repmat(repelem(bvec,na1,1),1,nz*np*nf);   % norefi: b'=b (same debt level)

b1rvec    = (1-theta)*exp(pvec); % refi: b'=(1-theta)p (at collateral constraint)
b1rinterp = repmat(repelem(b1rvec,na1,1),1,nz*np*nf); 
  
z1interp  = repmat(zmat',ns*na1,np*nf);
p1interp  = repmat(repelem(pmat',1,nz),ns*na1,nf);
f1interp  = repelem(fmat',ns*na1,np*nz);

zpf_prob   = repelem(kron(fprob,kron(ptran,ztran)),nb*na*na1,1);

% consumption and u
wlthvec = exp(zvec)+(1+ra)*avec-(1+rb)*bvec;  

cr = repmat(wlthvec+(1-fb)*b1rvec,1,na1)-a1choice;  %ns*na1
ur = u(cr,s,infi,mar) + repmat(fvec,1,na1);

cn = repmat(wlthvec+bvec,1,na1)-a1choice;    %ns*na1
un = u(cn,s,infi,mar);


% ============== value function solution  =============

Vold = Vini;
diff = diffini;
iter = 0;

while (diff>difftol)

    % Interpolate
    if fprob==1   % frictionless case
v1nvec = interpn(bmat,amat,zmat,pmat,Vold,b1ninterp,a1interp,z1interp,p1interp,'linear');
v1rvec = interpn(bmat,amat,zmat,pmat,Vold,b1rinterp,a1interp,z1interp,p1interp,'linear');
    else
v1nvec = interpn(bmat,amat,zmat,pmat,fmat,Vold,b1ninterp,a1interp,z1interp,p1interp,f1interp,'linear');
v1rvec = interpn(bmat,amat,zmat,pmat,fmat,Vold,b1rinterp,a1interp,z1interp,p1interp,f1interp,'linear');
    end
v1nvec(isnan(v1nvec)) = -infi;
v1rvec(isnan(v1rvec)) = -infi;

    % expectation
ev1nvec = sum(v1nvec.*zpf_prob,2);
ev1rvec = sum(v1rvec.*zpf_prob,2);

ev1n    = (reshape(ev1nvec,na1,ns))';
ev1r    = (reshape(ev1rvec,na1,ns))';

    % value function of all states, all choices
vn = un+beta*ev1n;  % ns*na1
vr = ur+beta*ev1r;  % ns*na1

[maxvr,maxvridx] = max(vr,[],2);
[maxvn,maxvnidx] = max(vn,[],2);

Vvec = max([maxvn maxvr],[],2);  % ns*1

Voldvec = reshape(Vold,ns,1);

diff = sum(abs(Vvec-Voldvec));

Vold = reshape(Vvec,nb,na,nz,np,nf);

iter = iter+1;

disp(iter)
disp(diff)
end

V = Vold;



% policy functions (based on the last iteration)

%adjvec = (maxvr>maxvn);
%ridx = sub2ind([ns na1],(1:ns)',maxvridx);
%nidx = sub2ind([ns na1],(1:ns)',maxvnidx);

%a1vec  = adjvec.*a1choice(ridx)+(1-adjvec).*a1choice(nidx);
%b1vec  = adjvec.*b1rvec  + (1-adjvec).*bvec;
%cvec   = adjvec.*cr(ridx)+(1-adjvec).*cn(nidx);

% convert to matrix

%adj  = reshape(adjvec,nb,na,nz,np,nf);
%a1   = reshape(a1vec,nb,na,nz,np,nf);
%c    = reshape(cvec,nb,na,nz,np,nf);

%b1   = reshape(b1vec,nb,na,nz,np,nf);
%b    = reshape(bvec,nb,na,nz,np,nf);
%refi = adj.*(b1>b);






